library(dplyr)
library(stats)
library(sjPlot)
library(tidyverse)
library(estimatr)
library(ggplot2)
library(plm)
library(stringr)
library(caTools)
library(car)
library(quantmod)
library(MASS)
library(corrplot)
library(brant)

####Load and join data####

#Set working directory
setwd('/Users/celinascott-buechler/DFP/dac_survey/Replication/Analysis/National survey/Data')

#Read in cleaned data file outputted from script 1
dac_results <- read_csv('dac_cleaned_final.csv', col_names = TRUE)

#Median income by zipcode downloaded from: https://data.census.gov/table/ACSST5Y2022.S1903?q=Income%20by%20Zip%20code%20tabulation%20area
mean_income <- read_csv('ACS_2022_Mean_Income.csv', col_names = TRUE)

mean_income$Zip <- substr(mean_income$Zip, start = 7, stop = nchar(mean_income$Zip))

dac_results <- left_join(dac_results, mean_income, by = c("ZipCode" = "Zip"))

#Tract to zipcode conversion downloaded from: https://www.huduser.gov/portal/datasets/usps_crosswalk.html#data
tract_zip <- read_csv('ZIP_TRACT_122021.csv', col_names = TRUE )

#CDC EJI indicators downloaded from: https://www.atsdr.cdc.gov/placeandhealth/eji/index.html
cdc_eji <-read.csv('CDC_EJI.csv')

#DOE energy community ID downloaded from: https://arcgis.netl.doe.gov/portal/apps/experiencebuilder/experience/?id=a2ce47d4721a477a8701bd0e08495e1d
coal_comms <- read.csv('Coal_Closure_Energy_Communities_2023v2.csv')

dac_results$HouseholdUnionization <- gsub("(.*),.*", "\\1", dac_results$HouseholdUnionization)
dac_results$HouseholdUnionization <- gsub("(.*),.*", "\\1", dac_results$HouseholdUnionization)
dac_results$HouseholdUnionization <- gsub("(.*),.*", "\\1", dac_results$HouseholdUnionization)
dac_results$HouseholdUnionization <- gsub("(.*),.*", "\\1", dac_results$HouseholdUnionization)

zipcodes <- read_csv('RUCA2010zipcode.csv', col_names=TRUE)

dac_results$ZipCode[dac_results$ZipCode == 20721-4503] <- 20721

dac_results <- left_join(dac_results, zipcodes, by = c("ZipCode" = "ZIP_CODE"))
dac_results$LocType[dac_results$RUCA1 == 1 | dac_results$RUCA1 == 2 | dac_results$RUCA1 == 3] <- 'Metropolitan'
dac_results$LocType[dac_results$RUCA1 == 4 | dac_results$RUCA1 == 5 | dac_results$RUCA1 == 6] <- 'Micropolitan'
dac_results$LocType[dac_results$RUCA1 == 7 | dac_results$RUCA1 == 8 | dac_results$RUCA1 == 9] <- 'Small town'
dac_results$LocType[dac_results$RUCA1 == 10] <- 'Rural'

#Use World Economic Forum income classifications to group income brackets: https://www.weforum.org/agenda/2022/07/household-income-distribution-wealth-inequality-united-states/
dac_results$Income_dup<- dac_results$Income
dac_results$Income_dup[dac_results$Income_dup=='Under $25,000'] <- '$0 - $25,000'
dac_results$Income_dup[dac_results$Income_dup=='More than $200,000'] <- '$200,000 - $300,000'
dac_results <- dac_results %>% separate(Income_dup, " - ", into=c("Income_lower", "Income_upper"))
dac_results$Income_lower <- as.numeric(str_replace_all(dac_results$Income_lower, "[^[:alnum:]]", ""))
dac_results$Income_upper <- as.numeric(str_replace_all(dac_results$Income_upper, "[^[:alnum:]]", ""))

dac_results$HouseholdSize[dac_results$HouseholdSize == '8+'] <- 8
dac_results$HouseholdSize <- as.numeric(dac_results$HouseholdSize)

dac_results$Income_lower_norm <- dac_results$Income_lower/dac_results$HouseholdSize*2.6
dac_results$Income_upper_norm <- dac_results$Income_upper/dac_results$HouseholdSize*2.6

dac_results$Income_med_norm <- apply(dac_results[, c("Income_lower_norm", "Income_upper_norm")], 1, median)

dac_results$HouseholdIncomeBracket <- NA

dac_results$HouseholdIncomeBracket[dac_results$Income_med_norm < 52000] <- 'lower'
dac_results$HouseholdIncomeBracket[dac_results$Income_med_norm > 52000 & dac_results$Income_med_norm < 156000] <- 'middle'
dac_results$HouseholdIncomeBracket[dac_results$Income_med_norm > 156000] <- 'upper'

dac_results$nonwhite[dac_results$Race!='White'] <- 1
dac_results$nonwhite[dac_results$Race=='White'] <- 0

dac_results$race[dac_results$nonwhite==0] <- 'white'
dac_results$race[dac_results$nonwhite==1] <- 'nonwhite'

dac_results$IncomeBracket <- "not low-income"
dac_results$IncomeBracket[dac_results$HouseholdIncomeBracket == "lower"] <- "low-income"


####Create numerical variables for regressions & t tests####

dac_results$PostSupportUS_num[dac_results$PostConjointSupportInUS == "Don't know"] <- NA
dac_results$PostSupportUS_num[dac_results$PostConjointSupportInUS == "Strongly oppose"] <- -2
dac_results$PostSupportUS_num[dac_results$PostConjointSupportInUS == "Somewhat oppose"] <- -1
dac_results$PostSupportUS_num[dac_results$PostConjointSupportInUS == "Somewhat support"] <- 1
dac_results$PostSupportUS_num[dac_results$PostConjointSupportInUS == "Strongly support"] <- 2

dac_results$PostSupportCommunity_num[dac_results$PostConjointSupportInCommunity == "Don't know"] <- NA
dac_results$PostSupportCommunity_num[dac_results$PostConjointSupportInCommunity == "Strongly oppose"] <- -2
dac_results$PostSupportCommunity_num[dac_results$PostConjointSupportInCommunity == "Somewhat oppose"] <- -1
dac_results$PostSupportCommunity_num[dac_results$PostConjointSupportInCommunity == "Somewhat support"] <- 1
dac_results$PostSupportCommunity_num[dac_results$PostConjointSupportInCommunity == "Strongly support"] <- 2

dac_results$PostSupportNearCommunity_num[dac_results$`PostConjointSupportNear Community` == "Don't know"] <- NA
dac_results$PostSupportNearCommunity_num[dac_results$`PostConjointSupportNear Community` == "Strongly oppose"] <- -2
dac_results$PostSupportNearCommunity_num[dac_results$`PostConjointSupportNear Community` == "Somewhat oppose"] <- -1
dac_results$PostSupportNearCommunity_num[dac_results$`PostConjointSupportNear Community` == "Somewhat support"] <- 1
dac_results$PostSupportNearCommunity_num[dac_results$`PostConjointSupportNear Community` == "Strongly support"] <- 2

dac_results$PostCommunitySupportUS_num[dac_results$PostConjointCommunitySupportInUS == "Don't know"] <- NA
dac_results$PostCommunitySupportUS_num[dac_results$PostConjointCommunitySupportInUS == "Strongly oppose"] <- -2
dac_results$PostCommunitySupportUS_num[dac_results$PostConjointCommunitySupportInUS == "Somewhat oppose"] <- -1
dac_results$PostCommunitySupportUS_num[dac_results$PostConjointCommunitySupportInUS == "Somewhat support"] <- 1
dac_results$PostCommunitySupportUS_num[dac_results$PostConjointCommunitySupportInUS == "Strongly support"] <- 2

dac_results$PostCommunitySupportCommunity_num[dac_results$PostConjointCommunitySupportInCommunity == "Don't know"] <- NA
dac_results$PostCommunitySupportCommunity_num[dac_results$PostConjointCommunitySupportInCommunity == "Strongly oppose"] <- -2
dac_results$PostCommunitySupportCommunity_num[dac_results$PostConjointCommunitySupportInCommunity == "Somewhat oppose"] <- -1
dac_results$PostCommunitySupportCommunity_num[dac_results$PostConjointCommunitySupportInCommunity == "Somewhat support"] <- 1
dac_results$PostCommunitySupportCommunity_num[dac_results$PostConjointCommunitySupportInCommunity == "Strongly support"] <- 2

dac_results$PostCommunitySupportNearCommunity_num[dac_results$`PostConjointCommunitySupportNear Community` == "Don't know"] <- NA
dac_results$PostCommunitySupportNearCommunity_num[dac_results$`PostConjointCommunitySupportNear Community`== "Strongly oppose"] <- -2
dac_results$PostCommunitySupportNearCommunity_num[dac_results$`PostConjointCommunitySupportNear Community` == "Somewhat oppose"] <- -1
dac_results$PostCommunitySupportNearCommunity_num[dac_results$`PostConjointCommunitySupportNear Community` == "Somewhat support"] <- 1
dac_results$PostCommunitySupportNearCommunity_num[dac_results$`PostConjointCommunitySupportNear Community` == "Strongly support"] <- 2

dac_results$PreSupportUS_num[dac_results$PreConjointSupportInUS == "Don't know"] <- NA
dac_results$PreSupportUS_num[dac_results$PreConjointSupportInUS == "Strongly oppose"] <- -2
dac_results$PreSupportUS_num[dac_results$PreConjointSupportInUS == "Somewhat oppose"] <- -1
dac_results$PreSupportUS_num[dac_results$PreConjointSupportInUS == "Somewhat support"] <- 1
dac_results$PreSupportUS_num[dac_results$PreConjointSupportInUS == "Strongly support"] <- 2

dac_results$PreSupportCommunity_num[dac_results$PreConjointSupportInCommunity == "Don't know"] <- NA
dac_results$PreSupportCommunity_num[dac_results$PreConjointSupportInCommunity == "Strongly oppose"] <- -2
dac_results$PreSupportCommunity_num[dac_results$PreConjointSupportInCommunity == "Somewhat oppose"] <- -1
dac_results$PreSupportCommunity_num[dac_results$PreConjointSupportInCommunity == "Somewhat support"] <- 1
dac_results$PreSupportCommunity_num[dac_results$PreConjointSupportInCommunity == "Strongly support"] <- 2

dac_results$PreSupportNearCommunity_num[dac_results$`PreConjointCommunitySupportNear Community` == "Don't know"] <- NA
dac_results$PreSupportNearCommunity_num[dac_results$`PreConjointCommunitySupportNear Community` == "Strongly oppose"] <- -2
dac_results$PreSupportNearCommunity_num[dac_results$`PreConjointCommunitySupportNear Community` == "Somewhat oppose"] <- -1
dac_results$PreSupportNearCommunity_num[dac_results$`PreConjointCommunitySupportNear Community` == "Somewhat support"] <- 1
dac_results$PreSupportNearCommunity_num[dac_results$`PreConjointCommunitySupportNear Community` == "Strongly support"] <- 2

dac_results$JobsNeed_num[dac_results$JobsNeed == "Don't know"] <- NA
dac_results$JobsNeed_num[dac_results$JobsNeed == "Little to no need"] <- 0
dac_results$JobsNeed_num[dac_results$JobsNeed == "Some need"] <- 1
dac_results$JobsNeed_num[dac_results$JobsNeed == "Large need"] <- 2

dac_results$Ideology_num[dac_results$Ideology == "Moderate"] <- 0
dac_results$Ideology_num[dac_results$Ideology == "Apolitical"] <- NA
dac_results$Ideology_num[dac_results$Ideology == "Something else"] <- NA
dac_results$Ideology_num[dac_results$Ideology == "Somewhat conservative"] <- -1
dac_results$Ideology_num[dac_results$Ideology == "Very conservative"] <- -2
dac_results$Ideology_num[dac_results$Ideology == "Somewhat liberal"] <- 1
dac_results$Ideology_num[dac_results$Ideology == "Very liberal"] <- 2

dac_results$HouseholdIncomeBracket_num[dac_results$HouseholdIncomeBracket == "lower"]<- 1
dac_results$HouseholdIncomeBracket_num[dac_results$HouseholdIncomeBracket == "middle"]<- 2
dac_results$HouseholdIncomeBracket_num[dac_results$HouseholdIncomeBracket == "upper"]<- 3

dac_results$ExtremeWeather_num[dac_results$ExtremeWeather == "Not at all concerned"] <- 0
dac_results$ExtremeWeather_num[dac_results$ExtremeWeather == "Not very concerned"] <- 1
dac_results$ExtremeWeather_num[dac_results$ExtremeWeather == "Somewhat concerned"] <- 2
dac_results$ExtremeWeather_num[dac_results$ExtremeWeather == "Very concerned"] <- 3
dac_results$ExtremeWeather_num[dac_results$ExtremeWeather == "Extremely concerning"] <- 4

dac_results$ClimateChange_num[dac_results$ClimateChange == "Not at all concerned"] <- 0
dac_results$ClimateChange_num[dac_results$ClimateChange == "Not very concerned"] <- 1
dac_results$ClimateChange_num[dac_results$ClimateChange == "Somewhat concerned"] <- 2
dac_results$ClimateChange_num[dac_results$ClimateChange == "Very concerned"] <- 3
dac_results$ClimateChange_num[dac_results$ClimateChange == "Extremely concerning"] <- 4

dac_results$AirWaterPollution_num[dac_results$AirWaterPollution == "Not at all concerned"] <- 0
dac_results$AirWaterPollution_num[dac_results$AirWaterPollution == "Not very concerned"] <- 1
dac_results$AirWaterPollution_num[dac_results$AirWaterPollution == "Somewhat concerned"] <- 2
dac_results$AirWaterPollution_num[dac_results$AirWaterPollution == "Very concerned"] <- 3
dac_results$AirWaterPollution_num[dac_results$AirWaterPollution == "Extremely concerning"] <- 4

dac_results$IndustryPromises_num[dac_results$IndustryPromises=="I expect industries will deliver both the number and quality of jobs they promise."] <- 2
dac_results$IndustryPromises_num[dac_results$IndustryPromises=="I expect industries will deliver the number but not the quality of jobs they promise." | dac_results$IndustryPromises=="I expect industries will deliver the quality but not the number of jobs they promise."] <- 1
dac_results$IndustryPromises_num[dac_results$IndustryPromises=="I expect industries will deliver neither the number nor the quality of jobs they promise."] <- 0
dac_results$IndustryPromises_num[dac_results$IndustryPromises=="Don’t know"] <- NA

dac_results$Party_num[dac_results$Party == "Independent"] <- 0
dac_results$Party_num[dac_results$Party == "Democrat"] <- 1
dac_results$Party_num[dac_results$Party == "Republican"] <- -1
dac_results$Party_num[dac_results$Party == "Something else"] <- NA

dac_results$CommunityReliance_num[dac_results$CommunityReliance == "Don't know"] <- NA
dac_results$CommunityReliance_num[dac_results$CommunityReliance == "Strongly disagree"] <- -2
dac_results$CommunityReliance_num[dac_results$CommunityReliance == "Somewhat disagree"] <- -1
dac_results$CommunityReliance_num[dac_results$CommunityReliance == "Somewhat agree"] <- 1
dac_results$CommunityReliance_num[dac_results$CommunityReliance == "Strongly agree"] <- 2

dac_results$CommunityBelonging_num[dac_results$CommunityBelonging == "Don't know"] <- NA
dac_results$CommunityBelonging_num[dac_results$CommunityBelonging == "Strongly disagree"] <- -2
dac_results$CommunityBelonging_num[dac_results$CommunityBelonging == "Somewhat disagree"] <- -1
dac_results$CommunityBelonging_num[dac_results$CommunityBelonging == "Somewhat agree"] <- 1
dac_results$CommunityBelonging_num[dac_results$CommunityBelonging == "Strongly agree"] <- 2

dac_results$CommunityMoveOut_num[dac_results$CommunityMoveOut == "Don't know"] <- NA
dac_results$CommunityMoveOut_num[dac_results$CommunityMoveOut == "Strongly disagree"] <- -2
dac_results$CommunityMoveOut_num[dac_results$CommunityMoveOut == "Somewhat disagree"] <- -1
dac_results$CommunityMoveOut_num[dac_results$CommunityMoveOut == "Somewhat agree"] <- 1
dac_results$CommunityMoveOut_num[dac_results$CommunityMoveOut == "Strongly agree"] <- 2

dac_results$CommunityNaturalEnvironment_num[dac_results$CommunityNaturalEnvironment == "Don't know"] <- NA
dac_results$CommunityNaturalEnvironment_num[dac_results$CommunityNaturalEnvironment == "Strongly disagree"] <- -2
dac_results$CommunityNaturalEnvironment_num[dac_results$CommunityNaturalEnvironment == "Somewhat disagree"] <- -1
dac_results$CommunityNaturalEnvironment_num[dac_results$CommunityNaturalEnvironment == "Somewhat agree"] <- 1
dac_results$CommunityNaturalEnvironment_num[dac_results$CommunityNaturalEnvironment == "Strongly agree"] <- 2

dac_results$CaringForNature_num[dac_results$CaringForNature == "Don't know"] <- NA
dac_results$CaringForNature_num[dac_results$CaringForNature == "Strongly disagree"] <- -2
dac_results$CaringForNature_num[dac_results$CaringForNature == "Somewhat disagree"] <- -1
dac_results$CaringForNature_num[dac_results$CaringForNature == "Somewhat agree"] <- 1
dac_results$CaringForNature_num[dac_results$CaringForNature == "Strongly agree"] <- 2

dac_results$AnthroClimateChange_num[dac_results$AnthroClimateChange == "Don't know"] <- NA
dac_results$AnthroClimateChange_num[dac_results$AnthroClimateChange == "Strongly disagree"] <- -2
dac_results$AnthroClimateChange_num[dac_results$AnthroClimateChange == "Somewhat disagree"] <- -1
dac_results$AnthroClimateChange_num[dac_results$AnthroClimateChange == "Somewhat agree"] <- 1
dac_results$AnthroClimateChange_num[dac_results$AnthroClimateChange == "Strongly agree"] <- 2

dac_results$LocalGovTrust_num[dac_results$LocalGovTrust == "Don't know"] <- NA
dac_results$LocalGovTrust_num[dac_results$LocalGovTrust == "Highly distrust"] <- -2
dac_results$LocalGovTrust_num[dac_results$LocalGovTrust == "Somewhat distrust"] <- -1
dac_results$LocalGovTrust_num[dac_results$LocalGovTrust == "Somewhat trust"] <- 1
dac_results$LocalGovTrust_num[dac_results$LocalGovTrust == "Highly trust"] <- 2

dac_results$StateGovTrust_num[dac_results$StateGovTrust == "Don't know"] <- NA
dac_results$StateGovTrust_num[dac_results$StateGovTrust == "Highly distrust"] <- -2
dac_results$StateGovTrust_num[dac_results$StateGovTrust == "Somewhat distrust"] <- -1
dac_results$StateGovTrust_num[dac_results$StateGovTrust == "Somewhat trust"] <- 1
dac_results$StateGovTrust_num[dac_results$StateGovTrust == "Highly trust"] <- 2

dac_results$FedGovTrust_num[dac_results$FedGovTrust == "Don't know"] <- NA
dac_results$FedGovTrust_num[dac_results$FedGovTrust == "Highly distrust"] <- -2
dac_results$FedGovTrust_num[dac_results$FedGovTrust == "Somewhat distrust"] <- -1
dac_results$FedGovTrust_num[dac_results$FedGovTrust == "Somewhat trust"] <- 1
dac_results$FedGovTrust_num[dac_results$FedGovTrust == "Highly trust"] <- 2

dac_results$BaselineHeardSupportInCommunity_num[dac_results$BaselineHeardSupportInCommunity == "Don't know"] <- NA
dac_results$BaselineHeardSupportInCommunity_num[dac_results$BaselineHeardSupportInCommunity == "Strongly oppose"] <- -2
dac_results$BaselineHeardSupportInCommunity_num[dac_results$BaselineHeardSupportInCommunity == "Somewhat oppose"] <- -1
dac_results$BaselineHeardSupportInCommunity_num[dac_results$BaselineHeardSupportInCommunity == "Somewhat support"] <- 1
dac_results$BaselineHeardSupportInCommunity_num[dac_results$BaselineHeardSupportInCommunity == "Strongly support"] <- 2

dac_results$CommunityUnionization_num[dac_results$CommunityUnionization=="Don’t know"] <- NA
dac_results$CommunityUnionization_num[dac_results$CommunityUnionization=="No workers or hardly any workers in my community are unionized."] <- 0
dac_results$CommunityUnionization_num[dac_results$CommunityUnionization=="Some workers in my community are unionized."] <- 1
dac_results$CommunityUnionization_num[dac_results$CommunityUnionization=="Most workers in my community are unionized."] <- 2

#edits all entries in cdc_eji$GEOID that are 10 units long (rather than 11 like the others) to add a 0 to the beginning
cdc_eji$GEOID <- as.character(cdc_eji$GEOID)
cdc_eji$GEOID <- stringr::str_pad(cdc_eji$GEOID, 11, side="left", pad="0")
coal_comms$geoid_tract_2020 <- stringr::str_pad(coal_comms$geoid_tract_2020,11, side='left', pad='0')

coal_comms$geoid_county_2020 <- stringr::str_pad(coal_comms$geoid_county_2020, 5, side="left", pad="0")
cdc_eji$geoid_cty_2020 <- substr(cdc_eji$GEOID,1,5)

ej_zip <- left_join(cdc_eji, tract_zip, by = c('GEOID'='tract'))

coal_comms$energy_comm <- 1

ej_jt_zip <- left_join(ej_zip, coal_comms, by=c('GEOID' = 'geoid_tract_2020'))
ej_jt_zip["energy_comm"][is.na(ej_jt_zip["energy_comm"])] <- 0

ej_jt_zip_summ <- ej_jt_zip %>% group_by(zip) %>% summarize(mean_SPL_EJI = mean(SPL_EJI, na.rm=TRUE), max_SPL_EJI = max(SPL_EJI, na.rm=TRUE), mean_energy_comm = mean(energy_comm, na.rm=TRUE), max_energy_comm = max(energy_comm), sum_energy_comm = sum(energy_comm))

#Join EJ indicators with survey results
dac_results <- left_join(dac_results, ej_jt_zip_summ, by = c('ZipCode'='zip'))

#recalculate EJI index percentiles
dac_results <- mutate(dac_results, EJI_prcntl = mean_SPL_EJI/max(mean_SPL_EJI, na.rm=TRUE))

####Descriptive Statistics of Survey Sample####
table(dac_results$Gender)/length(dac_results$Gender)
table(dac_results$Race)/length(dac_results$Race)
table(dac_results$HouseholdIncomeBracket)/length(dac_results$HouseholdIncomeBracket)
table(dac_results$Education)/length(dac_results$Education)
table(dac_results$Party)/length(dac_results$Party)
table(dac_results$AgeGroup)/length(dac_results$AgeGroup)


####Knowledge of DAC prior to survey####

table(dac_results$BaselineHeardDAC)/length(dac_results$BaselineHeardDAC)

dac_results$BaselineHeardDAC_num[dac_results$BaselineHeardDAC == "Nothing at all"] <- 0
dac_results$BaselineHeardDAC_num[dac_results$BaselineHeardDAC == "Only a little"] <- 1
dac_results$BaselineHeardDAC_num[dac_results$BaselineHeardDAC == "Some"] <- 2
dac_results$BaselineHeardDAC_num[dac_results$BaselineHeardDAC == "A lot"] <- 3

tab_model(lm(BaselineHeardDAC_num~Age+Education_num+Party+Gender+nonwhite+HouseholdIncomeBracket_num, data = dac_results, weights=nationalweight))

pre_heard <- dac_results %>% filter(BaselineHeardDAC!="Nothing at all") #450 respondents
table(pre_heard$`PreConjointSupportNear Community`)/length(pre_heard$`PreConjointSupportNear Community`)
table(pre_heard$PreConjointSupportInCommunity)/length(pre_heard$PreConjointSupportInCommunity)
table(pre_heard$PreConjointSupportInUS)/length(pre_heard$PreConjointSupportInUS)
table(pre_heard$Party)/length(pre_heard$Party)

pre_heard$PrePostChangeUS <- pre_heard$PostSupportUS_num - pre_heard$PreSupportUS_num
pre_heard$PrePostChangeCommunity <- pre_heard$PostSupportCommunity_num - pre_heard$PreSupportCommunity_num

pre_heard_rep <- pre_heard %>% filter(Party=="Republican")
pre_heard_dem <- pre_heard %>% filter(Party=="Democrat")
pre_heard_ind <- pre_heard %>% filter(Party=="Independent")

tab_model(lm(PrePostChangeUS ~ Party, data=pre_heard))
tab_model(lm(PrePostChangeCommunity ~ Party, data=pre_heard))


####ORDINAL LOGISTIC REGRESSION####

#Set ordinal outcome variables
dac_results$PostConjointSupportInUS_factor <- factor(dac_results$PostConjointSupportInUS,
                                                     levels = c("Strongly oppose", "Somewhat oppose", "Somewhat support", "Strongly support"),
                                                     ordered = TRUE)

dac_results$PostConjointSupportInCommunity_factor <- factor(dac_results$PostConjointSupportInCommunity,
                                                            levels = c("Strongly oppose", "Somewhat oppose", "Somewhat support", "Strongly support"),
                                                            ordered = TRUE)

dac_results$PostConjointSupportNearCommunity_factor <- factor(dac_results$`PostConjointSupportNear Community`,
                                                            levels = c("Strongly oppose", "Somewhat oppose", "Somewhat support", "Strongly support"),
                                                            ordered = TRUE)

#Construct model for in-community support
in_comm<- polr(PostConjointSupportInCommunity_factor ~ JobsNeed_num+ IndustryPromises_num+CommunityReliance_num+AnthroClimateChange_num * Party + Age + Education_num +IncomeBracket + Gender +nonwhite, data = dac_results, Hess=TRUE, weights = nationalweight)

#View model results
tab_model(in_comm)
summary(in_comm)

#Check proportional odds assumption using Brant's test
brant(in_comm)

#Check VIF for multicollinearity
vif(in_comm)


#Construct model for near-community support
near_comm<- polr(PostConjointSupportNearCommunity_factor ~ JobsNeed_num+ IndustryPromises_num+CommunityReliance_num+AnthroClimateChange_num * Party + Age + Education_num +IncomeBracket + Gender +nonwhite, data = dac_results, Hess=TRUE, weights = nationalweight)

#View model results
tab_model(near_comm)
summary(near_comm)

#Check proportional odds assumption using Brant's test
brant(near_comm)

#Check VIF for multicollinearity
vif(near_comm)


#Construct model for in-US support
in_us<- polr(PostConjointSupportInUS_factor ~ JobsNeed_num+ IndustryPromises_num+CommunityReliance_num+AnthroClimateChange_num * Party + Age + Education_num +IncomeBracket + Gender +nonwhite, data = dac_results, Hess=TRUE, weights = nationalweight)

#View model results
tab_model(in_us)
summary(in_us)

#Check proportional odds assumption using Brant's test
brant(in_us)

#Check VIF for multicollinearity
vif(in_us)


###Environmental Justice and Just Transition indicators####

tab_model(lm(mean_SPL_EJI~Party+Age+Education_num+Gender+nonwhite+HouseholdIncomeBracket_num, data=dac_results))

in_comm_ej_jt<- polr(PostConjointSupportInCommunity_factor ~ mean_SPL_EJI+mean_energy_comm+JobsNeed_num+ IndustryPromises_num+CommunityReliance_num+AnthroClimateChange_num * Party + Age + Education_num +IncomeBracket + Gender +nonwhite, data=dac_results, Hess=TRUE, weights = nationalweight)

tab_model(in_comm_ej_jt)
summary(in_comm_ej_jt)
brant(in_comm_ej_jt)

near_comm_ej_jt<- polr(PostConjointSupportNearCommunity_factor ~ mean_SPL_EJI+mean_energy_comm+JobsNeed_num+ IndustryPromises_num+CommunityReliance_num+AnthroClimateChange_num * Party + Age + Education_num +IncomeBracket + Gender + nonwhite , data = dac_results, Hess=TRUE, weights = nationalweight)

tab_model(near_comm_ej_jt)
summary(near_comm_ej_jt)
brant(near_comm_ej_jt)

in_us_ej_jt<- polr(PostConjointSupportInUS_factor ~ mean_SPL_EJI+mean_energy_comm+JobsNeed_num+ IndustryPromises_num+CommunityReliance_num+AnthroClimateChange_num * Party + Age + Education_num +IncomeBracket + Gender +nonwhite, data = dac_results, Hess=TRUE, weights = nationalweight)

tab_model(in_us_ej_jt)
summary(in_us_ej_jt)
brant(in_us_ej_jt)


dac_energy_comms <- dac_results %>% filter(max_energy_comm==1)
dac_non_energy_comms <- dac_results %>% filter(max_energy_comm==0)

#Not significant
t.test(dac_energy_comms$PostSupportCommunity_num, dac_non_energy_comms$PostSupportCommunity_num,paired = FALSE, conf.level = 0.95)

dac_ej_comms <- dac_results %>% filter(EJI_prcntl>=0.75)
dac_non_ej_comms <- dac_results %>% filter(EJI_prcntl<0.75)

#Not significant
t.test(dac_ej_comms$PostSupportCommunity_num, dac_non_ej_comms$PostSupportCommunity_num,paired = FALSE, conf.level = 0.95)
